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ABSTRACT 

MG 2016+112 is a quadmply imaged lens system with two complete images A and B and 
a pair of merging partial images in region C as seen in the radio. The merging images are 
found to violate the expected mirror symmetry. This indicates an astrometric anomaly which 
could only be of gravitational origin and could arise due to substructure in the environment or 
line-of-sight of the lens galaxy. We present new high resolution multi-frequency VLBl obser- 
vations at 1.7, 5 and 8.4 GHz. Three new components are detected in the new VLBl imaging 
of both the lensed images A and B. The expected opposite parity of the lensed images A and B 
was confirmed due to the detection of non-collinear components. Furthermore, the observed 
properties of the newly detected components are inconsistent with the predictions of previous 
mass models. We present new scenarios for the background quasar which are consistent with 
the new observations. We also investigate the role of the satellite galaxy situated at the same 
redshift as the main lensing galaxy. Our new mass models demonstrate quantitatively that the 
satellite galaxy is the primary cause of the astrometric anomaly found in region C. The de- 
tected satellite is consistent with the abundance of subhaloes expected in the halo from cold 
dark matter (CDM) simulations. However, the fraction of the total halo mass in the satellite as 
computed from lens modeling is found to be higher than that predicted by CDM simulations. 

Key words: gravitational lensing - quasars: individual: MG 2016+112 - cosmology: obser- 
vations 



1 INTRODUCTION 

Cold dark matter (CDM) simulations of hierarchical struc- 
ture formation predict a large population of small mass sub- 
haloes within the extended dark matter parent halo es of 

galaxies (e.g.. [Nava rro. Frenk & White! 1 19961 : iMoore et alj|l999l : 
IPiemand, Kuhlen & Madau,2007 ). The total mass fraction of the 
parent halo that is made up in substructures is thought to be 5-10 
per cent and the mass function of the substructures is dN/dM oc 
M~^'^ (e.g., Diemand et al. 2007). Thus, a crucial test of the 
CDM paradigm is finding these substructures around galaxies like 
our own Milky Way. Although it appeared that the number of 
satellites found around our own galaxy, observed as dwarf galax- 
ies, was down by a factor of 10-100 compared to predictions 
(e.g., lKlvDinrtIl] ri999: Moore et al. 1999), in recent years the 
number of known Milky Way dwarf galaxies has increased (e.g.. 



IZucker et al.ll2006l: iBelokurov et"ai]|2006j|bl) . However, there still 
exists a disagreement in the number of observed sa tellites com - 
pared with the predictions from simulations ( Koposov et al.ll2008h . 
A possible solution to this 'missing satellites' problem is if part of 
the substructure popu lation is purely dark and devoid of a lu minous 
baryonic component teullock. Kravtsov & Weinberg|l200lh . 

Gravitational lensing offers a powerful method to detect CDM 
substructures in distant galaxies {z < 1) since the phenomenon 
is sensitive to matter that is both dark and luminous. The effect 
of small mass substructures within the lens galaxy is to change 
the magnifications and/or the positions of the lensed images 
from what is expec t ed fro m a globally smooth mass distribution 
(Bradac et al. 2002", '2004'; Chiba 2002; D alai & Koch anek 200i 
Metc alf & Zhao .2002; .Keeton ,2003,; .Keeton. Gaudi & Fetters 



2005:; iKochanek. Schneider & WambsganssI |2004| ; IChen et al 
l2007l ; lMcKean et alj|2007l) . 
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systemic bv lOalal & Kochanekl ( l2002h found that the mass fraction 
of the substructure required to satisfy the anomalous flux-densities 
of the lensed images was 0.6-7 per cent of the parent halo (90 per 
cent confidence level). At first glance, this appeared to be consis- 
tent with the expected level of C DM substr u cture from simula- 
tions (~5-10 per cent). However, iMao et al] ( |2004|) showed that 
in the inner part of dark matter haloes, that is, within a few kpc 
from the centre, which is the region probed by the multiple im- 
ages of the background source, the expected mass fraction of the 
CDM substructure is only < 0.5 per cent. Recent CDM simula- 
tions have predicted an even smaller mass fraction of < 0.3 per 
cent toiemand etaP 2007). 

The discrepancy between the inferred level of CDM sub- 
structure from gravitational lensing and dark matter simulations 
presents a potential problem for the CDM paradigm. This has 
prompted a number of studies into those systems with the worst 
known flux and astrometric anomalies in order to investigate 
their cause. Interestingly, two systems with anomalous lensed 
image flux-densities (CLASS B2045-H265 and MG 0414-^0534), 
have been found to harbour a small mass (dwarf) compan- 
ion galaxy to the lens, which when included in the lens mod- 
els, reproduce the observed flux-den sities of the lensed images 
jRos et al.ir2000l : iMcKean et alj|2007h The g ravitational lens sys- 
tem MG 2016-Hll j^l ^Lawrence et al.lll984r) . which is the focus 
of this paper, shows the most extreme example of the positions 
of the lensed images differing from the predictions of a globally 
smooth mass distribution. This astrometric anomaly has also been 
attributed to a satellite galaxy that is nearby the main lens ga laxy 
( iKochanek. Schneider & Wambsganssll2004i : IChen et al.ll2007h . In 
this paper, we present new high resolution radio imaging of the 
gravitational lens system MG 20I6-I-1 12 between 1.7 and 8.4 GHz. 
The purpose of these observations is to test mass models for the 
lensing galaxy and investigate the contribution, if any, of CDM sub- 
structure. 

Our paper is arranged as follows. We first present a review 
of the previous observations and mass models for MG 20I6-I-1 12 in 
Section|2l The new radio observations are presented in Section[3] In 
Section|4l the multi-frequency data are used to determine the radio 
spectral energy distributions of the lensed images. Mass models 
for the system are investigated in Section|5] We discuss our results 
in Section|6]and present our conclusions in Section|7] We use the 
following values of the cosmological parameters Sim = 0.3, JIa = 
0.7, h = 0.7. 



2 MG 2016-1-112 

MG 2016-1-112 ^Lawrence et al.lll984h was the first gravitational 
lens system to be discovered from the MIT-Green Bank 6-cm sur- 
vey (Bennett et al. 1986). The system show ed three compact radio 
components (A, B and C) in VLA imaging ^Lawrence et alj[T984l) 
which nearly form a right-angled triangle on the sky. The separa- 
tion between radio components A and B is ~ 3.4 arcsec and the 
flux-ratio is Sb/Sa ~ 1- The third radio component C, which 



^ Gravitational lens systems with a radio-loud lensed source are better for 
studying flux and astrometric anomalies because radio observations are i) 
not affected by dust extinction or micro-lensing and ii) provide high resolu- 
tion astrometry for the lensed images from very long baseline interferome- 
try. 

^ The nomenclature refers to B1950 coordinates. However, all of the ob- 
servations presented in this paper refer to J2000 coordinates. 



lies ~ 2 arcsec south-east of component B, dominates the radio 
emission from the system; the flux-ratio of components C and A is 
Sc / Sa ~ 3. Observations at a higher angular resolution (~50 mas 
beam size) with MERLIN at 5 GHz found components A and B 
to be compact, however, component C was resolved into essen- 
tially two further components, an extended CI and compact C2, 
separa ted by ~ 112 mas in an east- west direction dCarrett et al.l 
Il994ah . Very Long Baseline Interferometry (VLBI) observations 
at 1.7 GHz found extended structure in the components A and B 
jHeflin etaljl99lh whereas in the European VLBI Network (EVN) 
observations, component C was resolved into a bead of four com - 
ponents, identified as CI I, CI2, C13 and C2 jOarrett et al]|l996h . 
VLBI 5 GHz observations confirmed the presence of extended 
structure in radio components A, B and C dKoopmans et aljr2002l : 
hereafter referred to as K02). These data also showed that the com- 
ponents C12 and CI3 in region C had similar extended morpholo- 
gies, consistent with a pair of merging images. The outer pair, CI 1 
and C2, of region C also had a similar morphology, but their surface 
brightness' were different. 

High resolution imaging with the Hubble Space Telescope 
(HST) found the optical counterparts of radio components A and 
B to h ave strong point-s ource emission, consistent with a lensed 
quasar jFalco et al J 19991 . CfA-Arizona Space Telescope LEns Sur- 
vey (CASTLES)). Coincident with radio component C was an ex- 
tended gravitational arc, which had a redder colour than the quasar 
emission from A and B ^Lawrence. Ne ugebauer & Matthew j 
,1993.) . Optical spectroscopy showed A and B to be at redshift 
3.273 and hav e similar narrow emiss ion line spectra, typical of a 
type-2 quasar i Lawrence et al.|[l984l) . The extended gravitational 
arc was also found to be at a redshift of 3.273, however, the 
emission line ratios of the arc were different from those found 
in lensed images A and B jYamada et al.ll200ll) . Finally, the op- 
tical imaging also detected a massive elliptical galaxy between the 
three radio components. The lens, called galaxy D, is at a redshift 
1. 01 ( Schneider et al. 1986 i) and has a s tellar velocity dispersion 
of 328 kms^^ (Koopmans & Treij|2002h . The wide separation of 
components A an d B suggested that th e lens may be part of a larger 
structure. Indeed, ISoucailetal1j200lh spectroscopically confirmed 
an over-density of six galaxies at the same redshift as galaxy D al- 
though, subsequent Chandra observations failed to find evidence 
for a cluster in the form of extended X-ray em i ssion |Chartas et al.l 
l200lh . Moreover, IClowe. Trentham. & Tonrvl j200ll) identified a 
weak lensing mass peak (3 a) ~ 64 arcsec northwest of the sys- 
tem. 

The similar radio and optical morphologies and spectra of 
components A and B confirm that they are gravitationally lensed 
images of the same background source. However, the nature of 
component C is not so clear, since the radio spectrum is flatter than 
those of lensed images A and B, the radio morphology appears to be 
different and the optical emission line ratios are not consistent with 
those found in the lensed images. A number of mass models have 
been proposed for MG 2016-1-1 12, which have region C being either 
part of the lens plane (e.g., Lawrence et al. 1984: Lawrenc e et al.l 
1993 ) or part of the b ackground quasar (e.g.. lSchneider et al.ll986l : 
^Langston et alj|l9"9ll) or belonging to both background and fore- 



ground (e.g.. lNarasimha & Chitrell989l : lNair & Garretlll997h . The 

most recent mass model was presented by K02, which had the caus- 
tic passing through radio component 2 of the background source 
(i.e., A2 and B2 in the image plane; see fig. 2 of K02). In this sce- 
nario, the radio components I and 2 of the background radio source 
and the optical quasar are doubly imaged, and part of component 
2 is quadruply imaged, with the third and fourth images being ob- 
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Table 1. A summary of the MERLIN, VLBI and HSA observations 



Array/ 


Total obs. 


Integ. 


Scan length (min) 


Frequency 


time (h) 


time (s) 


Target 


Phase calibrator 


MERLIN 










L7GHz 


14 


8 


3.7 


1.7 


5 GHz 


14 


8 


3.7 


0.8 


VLBI 










1.7 GHz 


17 


2 


4.5 


2.4 


5 GHz 


17 


1 


4.5 


1.5 


HSA 










8.4 GHz 


7.5 


1 


3.5 


2.3 



Table 2. The positions, peak surface brightness and total flux densities of 
MG 2016+1 12 from the MERLIN observations at 1.7 GHz. The '-' impUes 
that the integrated flux density could not be determined. 



Comp. 


RA 


Dec 


^peak 


^total 




(mas) 


(mas) 


(mJy beam^^) 


(mJy) 


A 


0±0.5 


0±0.5 


34.6±0.3 


34.9±0.5 


B 


-3005.8±0.5 


-1503.9±0.5 


34.8±0.3 


35.6±0.5 


Cla 


-2045.4±0.9 


-3246.2±0.9 


13.4±0.3 




Clb 


-2093.8±1.4 


-3221.5±1.4 


42.9±0.3 


62.8±0.6 



served in region C as the merging pair CI 2 and CI 3. This picture 
also explained the extended gravitational arc seen in the optical at 
region C, since here we are observing the quasar host galaxy with- 
out any contamination from the active core. Furthermore, if region 
C in fact corresponds to two merging images, then we would ex- 
pect the outer radio components CI 1 and C2 to have similar prop- 
erties. However, C2 is found to have a higher flux density. Also, 
the position of C2 is much farther to the west than expected. A 
possible solution to this ast rometric and flux-ratio anomaly is a 
nearby satellite galaxy, Gl jKochanek. Schneider & Wambsgansa 
[2004; Chen et al. 2007h. which is at the same redshift as the lens 
galaxy D dKoopmans & Treijl2002h . 



Table 3. The positions, peak surface brightness and total flux densities of 
MG 2016+1 12 from the MERLIN observations at 5 GHz. 



Comp. 


RA 


Dec 


^peak 


StotaL 




(mas) 


(mas) 


(mJy beam~^) 


(mJy) 


A 


0.0±0.5 


O.OitO.5 


11.7±0.2 


12.1±0.3 


B 


-3003.2±0.5 


-1502.0±0.5 


13.1±0.2 


13.1±0.3 


Cla 


-2048.3±0.5 


-3230.5±0.5 


15.4±0.2 


20.0±0.4 


Clb 


-2090.8±0.5 


-3225.6±0.5 


13.1±0.2 


13.2±0.3 


C2 


-2171.5±1.5 


-3212.8±1.5 


3.1±0.2 


3.0±0.3 



3 RADIO OBSERVATIONS 

Radio observations of MG 2016+1 12 are presented in this section. 
These include simultaneous MERLIN and global VLBI observa- 
tions of MG 2016+112 at both 1.7 and 5 GHz, and observations 
with the High Sensitivity Array (HSA) at 8.4 GHz. 



lens system. A phase referenced map of MG 2016+112 was sub- 
sequently made using IMAGR. Using the phase-referenced map as 
an initial model, phase self-calibration was performed on the lens 
system. The resulting map was used as a model for a new iteration 
of calibration of the data and this process was repeated until there 
was no significant improvement in the produced map. 



3.1 MERLIN 

3.1.1 Observations and data reduction 

The MERLIN array was used to observe MG 2016+112 simulta- 
neously with the global VLBI array at 1.7 and 5 GHz. The obser- 
vations were carried out on 2002 February 25 at 1.7 GHz and on 
2001 November 17 at 5 GHz. Both of the experiments had the same 
observational setup. The data were taken in a single band (referred 
to as intermediate frequency, IF) of 15 MHz bandwidth and sub- 
divided into 15 channels. The flux calibrator 3C 286 and the point 
source calibrator B2134+004 were observed for a few minutes each 
(one— two scans). The observations nodded between the lens sys- 
tem MG 2016+112 and the phase calibrator B2029+121. Further 
observational details can be found in Table[T] 

Using the standard MERLIN data reduction program, the data 
were first corrected for any non-closing errors and the amplitude 
calibration was performed. The remainder of the data calibration 
and imaging was carried out using the Astronomical Image Pro- 
cessing Software (AIPS) provided by the National Radio Astro- 
nomical Observatory (NRAO). With the phase referencing tech- 
nique, the fringe fitting solutions determined for the phase calibra- 
tor B2029+121 were applied to calibrate the phase information on 
MG 2016+112. Furthermore, amplitude and phase self-calibration 
was performed on the phase calibrator and applied to the target 



3.1.2 Results 

Fig. [T] shows the MERLIN maps of the lensed images of 
MG 2016+1 12 at 1.7 and 5 GHz. The weighting of the uv-data can 
be chosen between natural and uniform to give either better reso- 
lution or sensitivity at the cost of the other. It is possible however, 
to optimise the weighting by setting the parameter ROBUST to 
in IMAGR. Gaussian model fitting of the components in the image 
plane was done with the task JMFIT. The results are presented in Ta- 
bles[2]and[3] Lensed images A and B are found to be compact, and 
were each fitted by a single Gaussian model component at both fre- 
quencies. The slightly extended component C was fitted with two 
Gaussian components at 1 .7 GHz and with three at 5 GHz. Note that 
the components of CI that were detected at 1.7 and 5 GHz may not 
be the same. These components were identified by comparing with 
the model fitting of the high resolution VLBI images (see Section 
13. 2t . The relative flux densities of images B, CI (Cla+Clb) and 
C2 with respect to image A for the MERL IN 5 GHz imaging are 
consistent with previous measurements by iGarrett e t al. (1994al). 
However, the absolute flux densities of the lensed images in the 
new observations are less by 23 per cent which might be due to er- 
rors in the amplitude calibration of either the previous or new data 
sets. 
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Figure 1. Left: The MERLIN image of MG 2016-H112 at 1.7 GHz shows three lensed images (A, B and C). Lensed images A and B are compact whereas 
region C is sUghtly extended in the east-west direction. The map is restored with a beam of size 0.12x0.12 arcsec^. Right: The higher resolution 5 GHz 
MERLIN image shows that lensed images A and B are still compact, but region C is resolved into a brighter and extended component C 1 , and a compact 
component C2 to its west. Furthermore, the squares indicate the positions of galaxies D and Gl found in the optical. The map is restored with a beam of size 
0.07x0.04 arcsec^ and a position angle of 31 deg. The root-mean-square (rms) noise in the maps are 0.27 mjy beam~^ at 1.7 GHz and 0.18 mJy beam~^ at 
5 GHz. The contour levels are at (-3, 3, 6, 12, 24, 48) X the rms noise in the respective maps. North is up and east is to the left. 



3.2 Global VLBI 

3.2.1 Obsen'ations and data reduction 

Earlier 1.7 GHz EVN observations ( lGarrettetal1ll996l) revealed 
fine structure in region C. With a high resolution spectral analy- 
sis, the predictions of the spectra of the components in region C 
from some of the mass models (e.g, Nair & Garrett 1997) could 
be tested. Therefore, high resolution global VLBI observations of 
MG 2016-1-1 12 were undertaken at 1.7 and 5 GHz on 2002 February 
25 and 2001 November 17, respectively. 

Since the lensed images of MG 2016-1-112 have a low flux 
density, phase referenced observations are vital in determining the 
phase corrections on all baselines. A strong source within 2 deg of 
the lens system (B2029-I-121 with a total flux density ~ 0.9 Jy), 
was used as the phase calibrator for these observations. The scan 
lengths for each observation of MG 2016-1-112 and the calibrator 
B2029-F121, the correlator integration time and the total time of the 
observations are listed in Table [T] The data were taken in four IPs 
at 1.7 GHz and in two IPs at 5 GHz. Each IP had a bandwidth of 
8 MHz and was further divided into 16 channels. 

The antennas used for the 1.7 GHz observations were Effels- 
berg (Eb), Jodrell Bank (Jb), Medicina (Mc), Onsala (On), Torun 
(Tr) and the 10 antennas of the VLB A. In addition, the phased VLA 
(Y) and the 70 m dishes of Robledo and Goldstone were included 
for better sensitivity. Arecibo (Ar) and Westerbork (Wb) were also 
used, however, these data were lost at the VLBA correlator. At 
5 GHz, all of the VLBA antennas and Eb, Jb, Mc, On, Tr, Ar, Y 
and Wb were used. 

The data were calibrated in the standard manner within AIPS. 
The data were also corrected for the change in the parallactic an- 
gle to account for the apparent change in the position angle as the 
source moves across the sky. Subsequently, the a priori amplitude 



calibration was performed using the system temperatures of each 
antenna. To determine the system temperature of the phased VLA 
during the observations, the amplitude level of the phase calibrator 
was set to the known values obtained from the list of calibrators at 
different frequencies. 

The data for the calibrator were fringe fitted on the longest 
scan to determine the phase, delay and rate solutions for a sin- 
gle IF initially, and then for all of the IPs combined, which were 
then applied to the whole data set. Since averaging the data across 
the bandwidth can result in loss of amplitude and phase informa- 
tion in the wide-field maps, the data channels were not averaged, 
in spite of the large number of visibilities. The phase calibrator 
was mapped by performing phase self-calibration and then both 
phase and amplitude calibration. The phase and amplitude correc- 
tions from the best calibrated map of B2029-I-121 were then applied 
to the target. Subsequently, the data on the lens system were phase 
self-calibrated with solution intervals of 3 and 2 min for the 1.7 and 
5 GHz datasets, respectively. 

For MG 2016-1-112, the lensed sources within the field have 
a large angular separation (e.g., a few arcseconds). Hence, multi- 
ple small fields centred on the lensed images were mapped, that is, 
three windows centred at the position of A, B and region C were 
CLEANed while mapping MG 2016-1-112 with IMAGR. At 1.7 GHz, 
the weighting was chosen to produce an optimal balance between 
resolution and sensitivity, whereas at 5 GHz the resolution was 
good enough, so only the sensitivity had to be up-weighted, hence, 
natural weighting was applied. Similar to the technique used for 
the reduction of MERLIN data, self-calibration and mapping were 
iteratively carried out until no significant improvement was found 
in the maps. 
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3.2.2 The lensed images A and B 

Fig. E] shows the lensed images A and B, and the pair of merging 
images in region C at 1.7 GHz. At high resolution, the lensed image 
A shows a rich structure of five components. Image B also shows 
extended emission b ut is not so well resolv ed. In the earlier EVN 
5 GHz observations ('Koopma ns et al.ir200 2'). the lensed images A 
and B showed two components each with a hint of substructure. 
Since the pair of components AI-BI and A2-B2 have a high sur- 
face brightness, these were detected in the earlier EVN 5 GHz ob- 
servations. 

In the new global VLBI 1.7 GHz data presented here, images 
A and B are found to have three and two new components, respec- 
tively. Elliptical Gaussian model fitting of the images wa s carried 
out using Powell's minimization routine jPress et al]l 19921) and the 
results are given in Table|4] The errors at all of the frequencies were 
determined based on the principles described in Fomalont ( 1999). 
All of the components are numbered in decreasing order of their to- 
tal intensity. Image A is fitted with five Gaussian components and 
image B is fitted with only four Gaussian components. Since image 
B is the counterpart of image A with a relative magnification ~ 1, 
image B is also expected to have five components. The component 
B3-I-B2 is identified as the composite of the unresolved components 
B3 and B2. Furthermore, if B3 has a higher peak flux density than 
B2 (like its counterpart A3 in image A), then the peak position of 
B3-I-B2 will be closer to the peak position of B3. Thus, the model- 



fitted components of image A may not represent the counterpart 
components of image B. 

Fig. [3] shows the high resolution global VLBI images at 
5 GHz. Image A clearly shows the five components which were 
not well-resolved at 1.7 GHz. Here, image B also shows the ex- 
pected five-component structure, that is, with components B3 and 
B2 now clearly resolved. Component 2 in images A and B has a 
higher surface brightness than component 3, whereas at 1.7 GHz, 
component 3 has a higher surface brightness than component 2 in 
image A. This suggests the presence of frequency dependent struc- 
ture. The model fitting was again carried out with Powell's mini- 
mization routine. Images A and B were fitted with five components 
each, the results of which are presented in Table[5] 

The series of components observed in the lensed images A and 
B at 1.7 and 5 GHz are non-collinear which can be used to test the 
opposite parity expected from the theory of gravitational tensing. 
The opposite parity of images A and B is clearly demonstrated by 
the opposite curvature of the series of components. 

3.2.3 Region C 

Region C was previously known to have four components and all of 
these were detected in the 1.7 and 5 GHz imaging (see Fig. |2] and 
O. The pair of merging images in region C are referred to as CI 
(the east pair) and C2 (the west pair). The components have been 
numbered in an ascending order going inwards from the outside. 
The outer components are labelled as CI I and C21 for the east and 
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Figure 3. The global VLBI maps of lensed images A and B, and region C at 5 GHz with natural weighting. The rms noise in the maps is 0.06 mJy beam~^ 
and the contours are (—4, 4, 8, 16, 32) X the rms noise in the map. The size of the restoring beam is 3.7x 1.2 mas^ and the position angle is —7.54 deg. The 
grey-scales are in mJy beam"^ . North is up and east is left. 



west pair, respectively. The elongated components are labelled as 
C12 and C22, which are further resolved into several components. 
Each of these components are labelled as a, b etc. going inwards 
from the outside, for example, C12a, C12b, etc. in CI and C22a, 
C22b, etc. in C2 (see Fig.[3j- Note that this labelling convention is 
different from that used previously. 

At 1.7 GHz, the outer components on either side (Cll and 
C21) are fitted with two Gaussian model components each. How- 
ever, due to the low resolution and low signal-to-noise ratio (SNR), 
these are identified as a single component each, and their centroid 
(flux-density weighted) positions are reported here. The inner pair 
of elongated components were each fitted by three Gaussian model 
components. The peak positions, the peak flux densities and the to- 
tal intensities from the Gaussian model fitting are given in Table 
|4] Note that the components a, b and c found in component CI 2 
may not be the counterparts of C22 for three reasons. Firstly, these 
components are not well-resolved. Secondly, if these are a pair of 
merging images straddling the critical curve, then the magnification 
gradient is rapidly changing on either side of the critical curve. Ide- 
ally, it should change similarly on either side. However, practically, 
this may not be observed. Thirdly, the ratio of component separa- 
tions of the west pair C2 (C21-C22) and east pair CI (C11-C12) is 
asymmetric indicating different stretching on either side. 

At 5 GHz, the four components of region C are further re- 
solved and show fine structure. Therefore, centroid positions were 
determined for most of the components by fitting multiple Gaus- 
sian components. Here, component Cll is fitted with one Gaussian. 



C12 was modeled with six Gaussian components whereas C22 was 
fitted with seven Gaussian components. The component C21 is ex- 
tended, and hence, fitted with three unresolved components. Three 
centroid positions were determined for each of the elongated com- 
ponents C12 and C22, and one centroid position for C21. The peak 
positions, flux densities and total intensities for A and B are pre- 
sented in Table |5] The centroid positions, peak flux densities and 
the total flux densities are also given for region C in Table|5] 



3.3 High Sensitivity Array (HSA) 

3.3.1 Observations and data reduction 

Further high resolution observations at 8.4 GHz were needed to 
independently confirm the series of components detected in the 
lensed images at 5 GHz, and to better determine the spectra of 
lensed images A and B. Since the background quasar has a steeply 
falling radio spectrum (ai.7 = —0.94, where 5*1, oc u"), observa- 
tions at a higher frequency demanded increased sensitivity. More- 
over, to carry out a spectral analysis of the finely resolved structure, 
high frequency and high resolution imaging was needed. Therefore, 
MG 2016-^112 was observed with the HSA at 8.4 GHz. The HSA 
included the following large antennas: the 305m— Ar, 100m— Eb, 
100m— Green Bank Telescope and phased VLA, in addition to the 
10 VLB A antennas. 

The observations were made on 2006 April 30 and lasted for 
7.5 h. The right hand and left hand circular polarisation data were 
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Figure 4. The HSA maps of lensed images A and B, and region C at 8.4 GHz with a weighting between unifonn and natural. The rms noise in the maps is 
33 /iJy beam^^ and the contours are (—3, 3, 6, 12, 24, 48) X the rms noise in the map. The size of the restoring beam is 1.9x0.7 mas^ and the position angle 
is —7 deg. The grey-scales are in mjy beam" ^ . North is up and east is left. 



Table 4. The positions relative to component Al, peak surface brightness and total flux densities of the components of lensed images A and B, and region C 
at 1.7 GHz from the global VLBI observations. The phase-referenced position of component A 1 is 20'' 19"* 18. 188" + 11°27'14.638" for the epoch J2000 
at 1.7 GHz. 



Comp. 


RA 


Dec 


^peak 


^total 




(mas) 


(mas) 


(mJy beam~l) 


(mJy) 


Al 


().00±0.02 


().()0±0.02 


19.7±1.6 


24.3±2.5 


A2 


-12.2±0.1 


4.5±0.1 


2.4±0.5 


2.8±0.7 


A3 


-8.6±0.1 


3.8±0.1 


3.5±0.6 


4.5±0.9 


A4 


6.2±0.1 


-3.8±0.1 


4.1 ±0.6 


6.5±1.2 


A5 


-18.0±0.3 


5.5±0.3 


1.1±0.3 


1.7±0.6 


Bl 


-3005.74±0.02 


-1503.63±0.02 


11.5±0.1 


10.9±1.6 


B3-HB2 


-3011.0±0.1 


-1498.0±0.1 


3.1±0.6 


2.9±0.8 


B4 


-3004.5±0.1 


-1502.6±0.1 


8.5±1.0 


21.9±2.8 


B5 


-3012.5±0.2 


-1495.9±0.2 


2.5±0.5 


4.2±1.0 


cll 


-2013.6±0.2 


-3233.3±0.2 


2.1±0.4 


11.6±2.3 


cl2a 


-2045.4±0.1 


-3229.9±0.1 


4.1±0.6 


15.8±2.4 


cl2b 


-2053.2±0.3 


-3228.7±0.3 


1.5±0.4 


4.6±1.2 


cl2c 


-2061.0±0.3 


-3228.2±0.3 


1.5±0.4 


5.7±1.4 


c21 


-2178.8±0.1 


-3209.6±0.1 


4.9±0.7 


18.4±2.1 


c22a 


-2098.3±0.2 


-3221.3±0.2 


1.7±0.4 


7.2±1.7 


c22b 


-2092.2±0.3 


-3223.9±0.3 


3.1±0.5 


17.5±3.0 


c22c 


-2083.5±0.6 


-3227.3±0.6 


0.8±0.3 


3.3±1.2 
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Table 5. The positions relative to component Al, peak surface brightness and total flux densities of the components of lensed images A and B, and region C at 
5 GHz from the global VLBI observations. For all of region C, except Cll, the centroid positions and the total flux densities are given. The phase-referenced 
position of component Al is 20'' 19'" 18. 187" + 11°27'14.635" for the epoch J2000 at 5 GHz. 



Comp. 


RA 


Dec 


7 

^peak 


c 

'-'total 




(mas) 


(mas) 


(mjy beam ^ ) 


^'mT^/^ 

(mjyj 


Al 


0.00±0.02 


0.00±0.02 


4.4±0.6 


6.5±1.0 


A2 


-12.20±0.02 


4.44±0.02 


3.3±0.5 


3.6±0.7 


A'? 


— O. J^U. i 






1 s-t-rt ^ 
i .omu.D 


A A 


D.OitU. i 


— J. /itu. 1 


U.OltU.Z 


u.yitu.4 


A S 




7-i-n 
J. / mu. J 


n A-\-C\ 7 


i .omu. / 


Bl 


-3005.95±0.02 


-1503.94±0.02 


4.9±0.6 


7.1±1.1 


B2 


-3012.48±0.02 


-1497.20±0.02 


3.2±0.5 


3.0±0.6 


B3 


-3011.5±0.1 


-1499.4±0.1 


1.0±0.3 


2.2±0.6 


B4 


-3002.0±0.1 


-1505.7±0.2 


0.7±0.2 


1.5±0.5 


B5 


-3015.1±0.2 


-1492.7±0.2 


0.4±0.2 


O.7±0.3 


cll 


-2012.3±0.2 


-3234.2±0.2 


0.8±0.2 


4.7±1.4 


cl2a 


-2043.8±1.0 


-3231.0±1.0 


1.5±0.3 


4.8±1.1 


cl2b 


-2049.1 ±1.0 


-3230.3±1.0 


1.7±0.3 


9.1±0.3 


cl2c 


-2056.3±1.0 


-3229.6±1.0 


0.7±0.2 


4.6±1.1 


c21 


-2179.4±0.1 


-3210.4±0.1 


0.9±0.2 


5.2±1.0 


c22a 


-2097.6±1.0 


-3210.4±1.0 


1.3±0.3 


8.0±1.4 


c22b 


-2091.3±1.0 


-3224.9±1.0 


1.5±0.3 


7.3±1.3 


c22c 


-2086.2±1.0 


-3225.5±1.0 


0.9±0.2 


5.0±0.9 



recorded together in four IFs, each with 8 MHz bandwidth and 16 
channels. No cross polarisation was performed. The data were cor- 
related with the VLBA correlator using an integration time of 1 s 
to reduce time averaged smearing (further observational details can 
be found in Table [TJ. Ar had a power failure and problems with 
software which allowed observations for only 1.5 h from the 3 h 
window available on-source. The data were reduced by following a 
similar technique to that used in reducing the global 1.7 and 5 GHz 
data. The images were weighted to obtain an optimum combination 
of sensitivity and resolution in the maps as described in Sect. |3.1| 

3.3.2 Results 

The images of A, B and region C are shown in Fig. |4] The HSA 
imaging at 8.4 GHz confirmed the detection of three new com- 
ponents in images A and B, in addition to the two components 
known from the earlier observations. The model fitting of the im- 
ages A and B was done using the AIPS task JMFIT. A and B were 
both fitted with 5 Gaussian components each. The lensed images of 
MG 2016-1-112 are weak at 8.4 GHz (mJy level). Hence, the faint 
components 4 and 5 in images A and B and the pair of merging im- 
ages in C have a low SNR. The total flux densities of these compo- 
nents were therefore measured by integrating the surface brightness 
within the 3 a confidence boundary (where a is the rms map noise). 
The positions of these components were found from the peak of the 
surface brightness. The peak surface brightness, the integrated flux 
densities and the peak positions of the components are presented 
in Table [6] The errors were determined using the same method as 
before. 



3.4 Emission from a fifth lensed image or tlie lensing galaxy 

To estimate an upper limit on the flux density of a possible fifth 
image located in the vicinity of the lensing galaxy D, a fourth sub- 



field was mapped along with the sub-fields centred on images A, 
B and C. The sub-field centred on the optical position of the pri- 
mary lens D at 1.7, 5 and 8.4 GHz was CLEANed using the task 
IMAGR. The sizes of the fields mapped were 0.51x0.51 arcsec at 
1.7 GHz, 0.3x0.3 arcsec at 5 GHz and a field of size 0.2x0.2 arc- 
sec at 8.4 GHz. The fields were naturally weighted in order to 
achieve maximum sensitivity. However, no radio emission from a 
fifth component or from the lensing galaxy was found. The flux 
density limits (5a-level) are 0.41, 0.18 and 0.10 mJy beam~^ at 
1.7, 5 and 8.4 GHz, respectively. From lens models with a small 
constan t-density core-radius, a hi ghly demagnified image is pre- 
dicted ('N arasimha^ Subramanian. & Chitrelll986l) . Thus, the fifth 
lensed image is not expected to be detected in the observations pre- 
sented here. 



4 RADIO SPECTRAL ENERGY DISTRIBUTIONS 

Using the data from the high resolution radio observations of the 
lensed images, it is now possible to carry out a spectral analy- 
sis of the lensed radio components. In Fig. |5] the integrated spec- 
tral energy distributions of the two lensed images A and B, and 
of region C between 1.7 and 8.4 GHz is shown. Images A and 
B have similar radio spectra and flux densities, as expected from 
gravitational lensing, and the spectra are consistent with previous 
multi-epoch and multi-frequency observations. However, there is 
evidence for lensed image A having a lower flux-density than im- 
age B at 8.4 GHz. This is probably due to image A being more 
resolved at this frequency compared to image B (see from Fig. |4] 
that image A is extended in an east-west direction, which has the 
best angular resolution). The radio spectrum of region C is found to 
be slightly flatter relative to images A and B. Furthermore, the flux 
density of region C is significantly higher, as expected for highly 
magnified images near the critical curve. 

Also shown in Fig. [5] are the spectra of the five detected ra- 
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Table 6. The positions relative to Al. peak surface brightness and total flux densities of the components of lensed images A and B, and region C at 8.4 GHz 
from the HSA observations. The phase-referenced position of component Al is 20'U9'"18.187= + 11°27'14.634" for the epoch J2000 at 8.4 GHz. 



Comp. 


RA 


Dec 


^peak 


Stotal 




(mas) 


(mas) 


(mJy beam" ^ ) 


(mJy) 


Al 


O.OitO.l 


O.OitO.l 


1.6±0.1 


2.6±0.1 


A2 


-12.3±0.1 


4.5±0.1 


0.9±0.1 


1.4±0.1 


A3 


-8.8±0.1 


3.9±0.1 


0.3±0.1 


0.5±0.1 


A4 


6.4±0.4 


-3.7±0.4 


0.2±0.1 


0.3±0.1 


A5 


-17.9±0.2 


5.3±0.2 


0.1 ±0.1 


0.1±0.1 


Bl 


-3005.5±0.1 


-1503.7±0.1 


1.7±0.1 


2.8±0.1 


B2 


-3012.1±0.1 


-1497.0±0.1 


1.4±0.1 


1.8±0.1 


B3 


-3010.6±0.4 


-1500.3±0.4 


0.4±0.1 


0.6±0.1 


B4 


-3001.5±0.2 


-1505.9±0.2 


0.1 ±0.1 


0.3±0.1 


B5 


-3014.6±0.2 


-1492.5±0.2 


0.2±0.1 


0.3±0.1 


cll 


-2012.3±0.3 


-3233.3±0.3 


0.2±0.1 


0.4±0.1 


cl2 


-2047.8±0.3 


-3229.6±0.3 


0.6±0.1 


7.5±0.7 


c21 


-2177.8±0.5 


-3209.9±0.5 


0.2±0.1 


1.2±0.3 


c22 


-2092.8±0.3 


-3223.7±0.3 


0.5±0.1 


8.5±1.0 






Frequency (GHz) 

Figure 5. (left) The radio spectra of lensed images A and B, and region C. (centre) The radio spectra of the components of lensed image A. (right) The radio 
spectra of the components in region C. The data are taken from the global VLBl and HSA observations between 1.7 and 8.4 GHz. 



dio components of image A. Components 1, 3 and 4 are found to 
have steep radio spectra between 1.7 and 8.4 GHz, which is con- 
sistent with optically thin, extended radio structure. Component 5 
appears to have a flatter radio spectrum between 1.7 and 5 GHz, 
before steepening at 8.4 GHz. However, the uncertainties in the 
flux-densities of component 5 are quite large, so the overall radio 
spectrum could be steeper. The only genuine flat-spectrum radio 
component is 2, which shows a turnover at 5 GHz. This optically 
thick radio component is therefore assumed to be the core of the 
radio source and coincident with the strong optical/infrared quasar 
emission found in the HST imaging. Thus, like the optical core, 
component 2 is doubly imaged. 

Finally, Fig.|5]also shows the radio spectra of the four compo- 
nents in region C between 1.7 and 8.46 GHz. The merging images 
are expected to belong to the same part of the background quasar 
and hence, they should have similar spectra. The components of the 
inner pair (C12-C22) and the components of the outer pair (Cll- 
C21) show similar spectra within the uncertainties. The inner pair 
of elongated components of C (C12-C22) have flatter spectra than 



the outer pair of components (C11-C21). Since the components of 
the inner pair also have higher flux densities, these dominate the 
spectrum of region C at low resolutions, for example, as found in 
the MERLIN imaging presented in Section IbTI Also, there is ev- 
idence for the C2 components having a higher flux-density than 
their CI counterparts, which also adds to the case for there being 
a gravitational perturbation to the mass distribution, for example, 
from the nearby satellite galaxy (see Section[5](- 



5 MASS MODELS 

In this section, various mass models for the lens potential of 
MG 20I6-I-I12 are tested with increasing complexity. Given the 
large number of observational constraints and the uncertainty of 
the region C components, different combinations of observational 
constraints are also tested. The aim is to find the simplest scenario 
for the mass distribution and the background source that will fit the 
positions of the five radio components observed in images A and 
B, the flux density ratio of (at least) component Al and Bl, and the 
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positions of the four components in region C. Firstly, the predic- 
tions of the K02 mass model are tested for its consistency with the 
new observations. Next, two-galaxy mass models are studied which 
take into account all of the observed constraints and the dwarf com- 
panion to the lens. Lastly, three-galaxy models are also tested and 
the predictions of various mass models are further discussed. 

5.1 The K02 mass model 

The multi-frequency high resolution observations were conducted 
with the aim of testing the K02 mass model. As a sanity check, 
the K02 mass model was first reconstructed from the K02 con- 
straints. Note, however, that the K02 constraints did not include 
the outermost pair of components in region C (Cll and C21). 
H ence, the astrometric an omaly issue was not addressed in K02. 
In lKoopmans et aLll2002l the source component 2 was shown to 
be situated on a caustic. The positions of Al-Bl as doubly im- 
aged and A2-B2-C12b-C22b as quadruply imaged components 
along with the flux density ratio (Sb/Sa) were used. The mass 
model consisted of the main galaxy D (singular isothermal ellip- 
soid, hereafter SIE), a mass distribution Ml (singular isothermal 
sphere, hereafter SIS) which contributes to the convergence com- 
ing fr om the environment as detected in the weak lensing analy- 
sis of Iciowe. Trentham. & Toiirvl ( 1200 ll) . and a mass distribution 
M2 (SIS) due to another ph ysically nearby over- density of galaxies 
found spectroscopically bv lSoucail et al.l l l200lf) . The mass model- 
ing analyses carri ed out here we re done using the publicly available 
code GRAVLENS jXeetonllloOll) . 

In light of the rich core-jet structure found from the high reso- 
lution observations presented in this paper, the reconstructed mass 
model of K02 was used to predict the positions and flux-densities 
of the different components in the lensed images. This was done 
using the 5 GHz global VLBI data from Table|5] The observed and 
predicted image positions for all of the components in lensed im- 
ages A and B, and in region C are shown in Fig. [6] Here, the newly 
found components 3, 4 and 5 of image A were used to predict the 
positions of their counterparts. We find from this model that the 
components 4, 1, 3 and part of 2 are doubly imaged whereas com- 
ponent 5 and part of component 2 (c, o and e) are quadruply im- 
aged. The labels c, o and e correspond to the merging pair of com- 
ponents C11-C21, C12a-C22a and C12b-C22b respectiveljB An 
interesting result of this model is that component 5 is predicted to 
have four images, which are referred to as A5, B5, Cl-5 and C2-5. 
Components Cl-5 and C2-5 have opposite parity and are predicted 
to be about 100 mas on either side of region C (see Fig. |6l. Their 
magnification relative to A5 (or BSflis larger by a factor of ~ 10. 
Since gravitational lensing conserves the surface brightness of the 
lensed images, the counterparts of A5 and B5 in region C are ex- 
pected to be 10 times larger in solid angle. 

The detection of component 5 in region C is thus a crucial 
test of the K02 model. However, there is no evidence of emission 
near the expected position of component 5 in region C for any of 
the datasets presented here. In order to test that our MERLIN and 
VLBI datasets had sufficient resolution and surface brightness sen- 
sitivity to detect a possible fifth component in region C, we carried 

^ The labels c, o and e have no significance and are merely used for la- 
belling convenience. 

* Since the flux density ratio of image A and B is ~ 1 and the same is 
almost true between C 1 and C2, the components of A and C2 will be taken 
as the representative components in the discussion about magnification for 
the sake of simplicity. 



out simulations where artificial components were added to the uv- 
dataset and new maps were created. The MERLIN 5 GHz dataset 
was found to be capable of detecting component 5 (see Fig.|7]l. The 
non-detection of component 5 in region C in the actual MERLIN 
5 GHz map (c.f. Fig.[T) indicates that component 5 is not quadruply 
imaged and the K02 scenario is either incorrect, or that the surface 
brightness of component 5 in region C is lower than what is pre- 
dicted, and the K02 scenario needs to be modified to take this into 
account. 

It is noted that irrespective of the newly found discrepancy of 
the component 5 predictions with the K02 mass model, the K02 
mass model was originally not complete since it did not fit the po- 
sitions of the outer components (Cll and C21) on either side of 
region C. Also, from Fig.|6] it is clear that the observed positions 
of the three new components found here (component 3, 4, and 5) 
are not reproduced in the lensed images A and B by the K02 mass 
model. 

5.2 Constraints and priors 

The constraints on the positions of the lensed images are taken from 
the high resolution global VLBI data at 5 GHz (see Table |5). The 
mass models will not be sens itive to an astrometric shift of the com- 
ponents which are ^ 1 mas jKochanek. Schneider & WambsganssI 
(2004). Thus, the astrometric uncertainties of all components are set 
to a minimum of 1 mas, in spite of the higher precision offered by 
the VLBI observations. For all mass models, the flux density ratio 
Si3i/SAi is constrained to 1.09±0.22 (Tabled. 

The lensed images A and B show opposite parity and the same 
five components. Therefore, all five components are at least doubly 
imaged. The components in region C appear to consist of two op- 
posite parity merging images with each showing a pair of compact 
and extended source components. The optical emission line spectra 
confirms that region C is at the same redshift as the lensed images 
A and B. Therefore, it is unlikely that the components of region 
C correspond to a source other than the lensed quasar (images A 
and B), which happen to lie close enough to give the expected four- 
image configuration of a single object. Hence, region C is almost 
certainly related to the same background source, but it is not clear 
to which part observed in image A (or B), or anything that is unseen 
in image A (or B). 

The main lensing galaxy D and the other line-of-sight galax- 
ies are detected in the optical and their positions, measured with 
respect to the compact optical emission from image A, are taken 
from the CASTLES database. As radio component 2 in the lensed 
images A and B shows the flattest radio spectrum of all the de- 
tected components, it is assumed to be the radio core and coincident 
with the optical quasar emission (Fig.[5j- The ellipticity of the lens 
galaxy D in the optical is 0.43±0.01 with a position angle of —59 
deg east of north. Since galaxy D is a giant elliptical with a stel- 
lar velocity dispersion of ~ 328 km s~^, it is expected to have the 
largest contribution to the image splitting and the ellipticity of the 
halo. Thus, we apply a prior on the ellipticity and fix the position 
angle from the values obtained from the optical surface brightness 
profile of galaxy D. This is to make sure that the minimiza- 
tion does not converge to any unreasonable mass models with an 
otherwise lower as compared to a desired model. Finally, from 
the combined lensing and stellar dynamics analysis conducted by 
iTreu & Koopmansl hmj) . the logarithmic slope of the density pro- 
file of the lensing galaxy was found to be 7' = 2.0±0. 1 . Therefore, 
only elliptical isothermal mass profiles for galaxy D were consid- 
ered. 
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Figure 6. A comparison of the observed positions (black) of the lensed components in images A and B, and region C with those predicted (magenta) from 
the K02 mass model. The critical curve in the image plane (red) is shown to split region C. The panel S on the left shows the source positions, relative to the 
tangential caustic (blue) whereas the same region at larger scales is shown in the panel to the right. The axes are in arcsec, the offsets are relative to component 
Al. North is up and east is left. 
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Figure 8. The various scenarios A, B, C and D (from left to right) for the critical curve splitting the background source. As the critical curve move left to right, 
the multiplicities of the components detected in the high resolution VLBI imaging will change. 
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Figure 7. MERLIN 5 GHz image with component 5 artificially added 

to the i<\'-data. The image has been restored with a beam of size 
0.07 x 0.05 arcsec^ and position angle of 31 deg. 



5.3 A two-galaxy model 

We first test a two galaxy model for the mass distribution that in- 
cludes the dominant elliptical galaxy D (SIE) and the nearby dwarf 
galaxy Gl (SIS), which is spectroscopically confirmed to be at the 
same redshift as galaxy D. The offsets of galaxies D and Gl from 
the radio component A2 (assumed to be coincident with the optical 
emission from image A, see Fig.[T]l are ARA = — 1740±3 mas, 
ADec = -1782±3 mas, and ARA = -2499±28 mas, ADec = 
— 4037±28 mas, respectively. Including the satellite Gl is essential 
because of its proximity to region C, which shows the asymmetry 
in the separated pair of opposite parity features. An external shear 
is also included in the mass model. 

We show different scenarios of the background source strad- 
dling the tangential caustic in the Fig. [8] which are investigated for 
the two-galaxy model. The four possible scenarios are. Scenario 
A - the caustic goes through source component 2, Scenario B - the 
caustic goes between source components 2 and 5 such that it grazes 
source component 5, Scenario C - the caustic goes through source 
component 5, and Scenario D - the caustic is situated beyond source 
component 5. The results of the mass modeling for each scenario 
are summarized in Table|7] 

5.3.1 Scenario A 

Constraints: In this scenario, the caustic goes through source com- 
ponent 2 as was also the case for the K02 scenaric[fl The source 
components 4, 1, 3 and a part of 2 fall in the doubly imaged region, 
whereas a part of component 2 and the whole of component 5 are 
quadruply imaged (see Fig.[8]l. The inner components of C (i.e., la- 
belled - o and e) are chosen as the counterparts of A2 and an unseen 
A2e to the west, respectively, while the outer components (i.e., la- 
belled - c) are associated with some unseen component A2c. The 
uncertainties for all observed components are chosen to be 1 mas, 

^ Note that scenario A and the scenario of K02 are the same. The difference 
lies in the data constraints that were used to make the mass models. Hence, 
they predict slightly different image positions and magnifications. 
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Figure 9. The critical curves (red) and the caustics (green) for the two- 
galaxy lens model. The position of the source (S) relative to the tangential 
caustic is shown. The x and y axes show separations relative to component 
Al which is chosen as the origin throughout. 



whereas for the unseen components they are chosen to be 5 mas 
so that they do not contribute significantly to the x^- Exceptions to 
this are the uncertainties of the counterparts of component A5 in 
region C, that is, Cl-5 and C2-5. Their uncertainties are chosen to 
be very high (i.e., 10* mas) since their positions might be affected 
due to substructure in an unexpected way. A total of 41 constraints 
exist for this mass model. There are 14 free source positions and 
four free model parameters, that is, the Einstein radii of the two 
galaxies, the external shear and its position angle. The degrees of 
freedom (dof) are 23. 

Results/Predictions: The best-fitting model has a reduced ~ 3.5. 
The largest contribution to the total comes from the image po- 
sitions (x^ ~ 60, most of which is from the components of re- 
gion C) and the second largest is from the galaxy positions (total 
X^ ~ 18; most of which is due to the satellite Gl). The model- 
predicted Einstein radii of the galaxies D and Gl were found to be 
1.570 and 0.143 arcsec, respectively. The recovered ellipticity of 
D was 0.42. The external shear was found to be 10 per cent with 
a position angle of —41.5 deg measured East of North. The criti- 
cal curves and the caustics for the two-galaxy model can be seen 
in Fig. |9] The observed and fitted image positions are shown in 
the top and bottom panels of Fig. [TO] The relative magnification 
of the components in region C, with respect to A2-A2o, is very 
high. For example, Se/SA2e ~ 10^ and SJSa2c ~ 10^. Hence, 
the counterparts of region C in images A and B are predicted to 
be 100-1000 times fainter (at a level of ~10 /iJy) and 10-30 times 
smaller in size. These counterparts can not be detected with the sen- 
sitivity of the observations undertaken at any of the three frequen- 
cies presented here. However, component 5 is predicted to have two 
counter-images in region C which are referred to as Cl-5 and C2-5. 
The predicted flux density ratio is Sc2-5/Sa5 ~ 15, and a size that 
is at least four times larger than that of A5. Since Cl-5 and C2-5 
are not detected in the MERLIN observations, this scenario is ruled 
out. 

5.3.2 Scenario B 

Constraints: In this scenario, although the caustic is situated be- 
tween component 2 and 5, it is closer to component 5 such that C 1 1- 
C21-A5-B5 are its four lensed images whereas the inner elongated 
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Figure 10. The observed (black) and predicted (magenta) positions of the components in lensed images A and B, and region C for scenario A. The critical 
curve is shown in red. The axes are in arcsec and the offsets relative to component Al. North is up and east is left. 



pair of components o and e are associated with unseen components 
A5o and A5e to the east of component A5. Thus, components 4, 
1, 3 and 2 will be doubly imaged as shown in the second panel of 
Fig. [8] The number of constraints are 41 and the free parameters 
remain the same, thereby, giving a dof of 23. 



Results/Predictions: Not surprisingly, the reduced is 3.6 because 
this scenario is a slight modification of scenario A and is not ex- 
pected to modify the global mass model significantly. Furthermore, 
the individual contributions and the best-fitting model param- 
eters are also similar. However, the predictions are expected to 
change here. The positions of the components in images A and B 
are fitted within 1 mas except for components 4 and 5, which are 
fitted within 1.5 mas. The observed and model predicted image po- 
sitions are shown in the bottom panel of Fig.[TT] The relative mag- 
nifications of the C components with respect to A are predicted to 
be Se/SASe ~ 10^ and SJSas ~ 200. Since the inner components 
of C (o and e) have very high magnification, their counterparts in 
images A and B would be unseen. This is consistent with the obser- 
vations. However, the counterparts of component c at the position 
of A5 (or B5) are predicted to be ~ 100 times fainter whereas the 
observed relative magnification is ^ 10. Due to this inconsistency, 
this scenario is not acceptable either. 



5.3.3 Scenario C 

Constraints: Here, the caustic goes through component 5 such that 
the inner elongated pair (o-e) of region C are associated with A5 
and an unseen component to the west (A5e) that are both quadruply 
imaged. The pair of components c are associated with an unseen 
component A5c that would lie further to the north-west. Compo- 
nents 4, 1, 3 and 2 are doubly imaged in this scenario. The total 
number of constraints and free parameters are the same as before, 
hence the dof are 23. 

Results/Predictions: The best-fitted model parameters are the same 
as before. The reduced x^ is 3.6 and the individual x^ contributions 
are similar to the previous scenarios. The doubly imaged compo- 
nents 1, 2 and 3 are fitted within the 1 mas uncertainties whereas 
component 4 is fitted within 1.5 mas. The counterparts of c, o and 
e in image A (or B) are predicted to be within 1.5 mas west of 
the peak position of component A5. The relative magnification of 
region C is predicted to be 100-1000 times higher than their coun- 
terparts in component A5 (or B5). Therefore, the corresponding 
components in images A and B would have flux densities of the 
order 5-50 ^Jy, which is much lower than the rms noise level of 
the global VLBI 5 GHz observations. Moreover, it is certainly not 
possible to resolve these components since they would have sizes 
10—30 times smaller. Nevertheless, since these predicted compo- 
nents could not be detected, their possible existence is not incom- 
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patible with our new observations. However, tliis scenario may or 
may not stand true with better observational constraints in the fu- 
ture. 



5.3.4 Scenario D 

Constraints: In this scenario, the caustic is situated to the west of 
component 5 such that components 4, 1, 3, 2 and 5 are all doubly 
imaged and the components of C are associated with some unde- 
tected components in images A and B. The total number of con- 
straints is 45. There are 16 free source positions and 4 free param- 
eters and hence, the dof is 25. 

Results/Predictions: The reduced for the best-fitting model is 
~ 3.4 and the model parameters are similar to those given in sce- 
nario A. The doubly imaged components 1, 2 and 3 are fitted within 
their uncertainties (1 mas) whereas components 4 and 5 are fitted 
within 1.5 mas. The c, o and e components are predicted to have 
four images. Their counterparts in A and B are expected to have a 
relative demagnification of 100-1000. Such components would not 
have been detected in the multi-frequency observations presented 
here. Therefore, this scenario is also consistent with the observa- 
tions. In fact, observations that have 1000 times better sensitivity 
are needed to fully test this scenario. 



5.3.5 No constraints on the position of the lens mass components 

The cluster associated with the main galaxy D is believed to be a 
proto-cluster which is not centrally concentrated yet. The conclu- 
sion that the cluster is not v irialized is inferred from the absence of 
any diffuse X-ray emission jChartas et alj200ll) . Therefore, the op- 
tical position of the BCG (galaxy D) may not be coincident with the 
centre of the cluster. Mass models with no constraint on the position 
of the lensing galaxy were hence tested for different scenarios. The 
following results are described for scenario D. Initially, the position 
of the main galaxy D was allowed to be free. The best-fitting model 
shifted galaxy D by ~ 60 mas to the west. The reduced of 1.5 
was mainly improved by better fitting the image positions. The total 
of the galaxy position did not change significantly, which arises 
here solely due to the satellite Gl. The best-fitting parameters are 
similar to those obtained from the mass model with the position of 
galaxy D constrained. The Einstein radii are 1.551 and 0.14 arcsec 
for galaxies D and Gl, respectively. The fitted ellipticity is 0.43 and 
the external shear is 11 per cent with a position angle of —41.8 deg. 

5.4 A three-galaxy model 

In spite of the high external shear (~ 10 per cent) for the two 
galaxy mass model (SIE+SIS), the image positions of all compo- 
nents could not be reproduced within their uncertainties. Given 
the dense environment around the lensing galaxy and a handful 
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of line-of-sight nearby galaxies, an additional galaxy, G2 (SIS) 
was added to the mass model (i.e., SIE+SIS+SIS). The position 
of G2, relative to radio component A2 (assumed to be coincident 
with the optical emission from A) is ARA — — 5749±10 mas and 
ADec — +1759±10 mas. This three-galaxy model was tested for 
various scenarios. However, introducing G2 to the model did not 
result in any significant change in either the reduced or the best- 
fitted model parameters (see Table|7]for details). 



6 DISCUSSION 

In this section, we discuss the role of the satellite galaxy and 
the abundance of substructure in MG 2016-1-1 12. Furthermore, we 
compare our results for the substructure with the expectations of 
standard ACDM cosmological simulations. 



6.1 The satellite galaxy Gl 

The mass distribution in MG 2016-1-1 12 definitely consists of mul- 
tiple components. With a smooth mass model of SIE-l-shear centred 
at D, it is possible to fit the positions of the components only in 
images A and B. Fitting the components of region C results in a 
very poor fit with a reduced of about 180. Including the satellite 
galaxy (Gl), significantly improves the fit to all of the astromet- 
ric constraints including the astrometric anomaly found in region C 
(reduced ~ 3.5). Such a model allows us to estimate the mass 
of the satellite galaxy. 

For the SIE-l-SIS-l-shear models, the Einstein radius of the 
satellite galaxy (i.e., the SIS component) is found to be ~ 0.14 arc- 
sec. Hence, the mass of the satellite (within its Einstein radius) is 
found to be 8 x 10^ Mq whereas the mass of the primary lens D 
(within its Einstein radius) is lO'^^ Mq (see Eas. lAll and lA2t . Thus, 
the model predicts the fraction of the mass in the satellite galaxy to 
be 0.8 per cent. This is con sistent with the satellite mass fraction 
(1 per cent) suggested by iKochanek. Schneider & Wambsgans^ 
( I2OO4I) . The model-predicted velocity dispersion of the satellite 
galaxy is 99 kms~^ . This is consistent with the mass modeling 
results of lChen et al] ( 12007 '), who find the velocity dispersion to be 
87 < cr < 101 kms"^ for the satellite Gl. 



6.2 Abundance of subhaloes 

In this section, we test the predictions of substructure from 
CDM simulations vis-a-vis the luminous substructure detected in 
MG 2016-1-112. First, we determine the number of subhaloes with 
a mass greater than or equal to that of the satellite galaxy expected 
from simulations. Second, we compare the substructure mass frac- 
tion at the projected separation of the satellite galaxy with that ex- 
pected from simulations. 

For the first test, we use the results of iGao et alj ( l2004h who 
studied the subhalo populations over a wide range in halo mass us- 
ing high resolution simulations in a ACDM universe. They found 
that the mass function of subhaloes per unit halo mass can be ex- 
pressed as 



diV 
dm 



lO-''(m,,b/i/M0)-'''hM-V 



(1) 



The mass modeling results yield the total halo mass of the lens 
galaxy within the virial radius as A/vir ~ 2.1 x 10^^ Mq (see 
Eo. lA4t and the mass of the satellite galaxy within its Einstein ra- 
dius as 8 X 10" Mq. The number of subhaloes with mass a greater 



than or equal to the mass of the satellite galaxy (obtained by inte- 
grating Eq.[T]above the satellite galaxy mass) is ~ 3. The detection 
of a single satellite galaxy suggests that the observations are consis- 
tent with the predictions within 2 a. Note that the observed subhalo 
is strictly a lower limit. 

For the second test, we use the results of iMao et alj ( l2004h . 
who investigate the substructure mass fraction, faub, in haloes at 
low redshift from ACDM simulation^. They estimate /sub as a 
function of increasing projected distance within annuli of equal log- 
arithmic width in projected separation. The satellite galaxy Gl in 
MG 2016-1-112 is at ~0.05 r^ii of the main lensing galaxy D. The 
simulations predict /sub ~ 0.01 ± 0.0 1 at this p rojected separation 
(see the bottom left panel of fig. 1 of Mao et al. 2004). While cal- 
culating the substructure mass, , Mao et al.. (2004) only use the mass 
within 0.025 r^ii of the subhalo centre. For a fair comparison, we 
calculate the mass of Gl within 0.025 rvir using the modeling re- 
sults. The mass in the smooth halo component within the requisite 
annuli is found by integrating the projected surface density (see 
Eqs. |A5l & |A6t of the main halo appropriately. Thus, the fraction 
calculated from lensing, /sub = 0.09 (see Eq. lA7l (, is higher than 
the fraction predicted by the simulations. Nevertheless, these values 
are again consistent within Scr. 

It should be understood that the comparison of the mass frac- 
tion from lensing with that of simulations is biased due to several 
reasons. The simulated subhaloes are dark matter only whereas the 
substructure in MG 2016-1-112 is luminous and hence, consists of 
baryons. Predictions from lensing will tend to be biased towards the 
high end of the subhalo mass function because massive subhaloes 
will produce easily dete ctable observable effects (also see e.g., 
iBrvan, Mao & Kavll2008l) . Furthermore, low statistics is severely 
affecting the current comparison and hence, more examples simi- 
lar to MG 2016-1-1 12 and MG 0414-1-0534 are needed that could be 
meaningfully compared to the predictions of simulations. 



6.3 CDM substructure near images A and B 

In MG 2016-1-1 12, the lensed images A and B show similar spectra. 
They exhibit the expected opposite parity. A local perturbation in 
any of the images can be found by comparing the relative magni- 
fication matrix with the flux ratios of the images. The presence of 
non-coUinear structure in the images allows an accurate determi- 
nation of the relative magnification matrix. The comparison shall 
further strengthen our belief in the absence of substructure close to 
images A and B. 

The relative magnification (fir) of the lensed images, that is, 
the magnification of an image (B) with respect to another image 
(A), is simply the ratio of individual magnifications /is and fiA- 
Ideally, the surface brightness of the lensed images is conserved, 
hence, the observed ratio of the flux densities of the images (de- 
noted by S) should be equal to the ratio of their solid angles (de- 
noted by Lj). This can be expressed as, 

Sb i^B I , ,/,^ASm 
/iA OA iJJA 



(2) 



where the relative magnification matrix (M^j ) is defined as, 



Bly 



Mil M12 
M21 M22 



Al:, 

Aly 



Ml^Ai.O) 



^ The results do not evolve significantly with redshift (Andrea Maccio, 
priv. comm.) . 
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Table 7. The results from the mass modeling analysis and the fitted parameters. The xfp ™d xff denote the due to the image positions and flux densities, 
respectively. The Xgp denotes the due to the positions of all of the galaxies in a model, whereas X^th denotes the contribution from other priors, for 
example, the ellipticity. See text for further details. 



Model-Scenario 


Xtot 










Model parameters 


Comments/Status 


SIE+SlS+shear 
















Scenario A 


81.0 


60.0 


1.9 


18.0 


0.3 


b£)=1.570, P£)=(- 1-746. -1-777). 


predicts extra components in region C 














60=0.42, 7s= (0.10, -41.5), 


- not found in data 














f)Gi=0.143 , Pgi=(-2.572, -3.972) 


(Not acceptable) 


Scenario B 


83.4 


60.6 


1.9 


20.8 


0.05 


bo=1.568, Po=(- 1.746, -1.776), 


relative magnification of component 5 














60=0.43, 7s= (0.10, -41.5), 


component C21-do not match data 














bGi=0.146, Pgi=(-2.575, -3.976) 


(Not acceptable) 


Scenario C 


83.5 


61.5 


1.9 


20.0 


0.07 


bi3=1.569, Po=(- 1.746, -1.776), 


relative magnifications are consistent 














60=0.43, 7s= (0.10, -41.5), 


with the data presented here 














f)Gi=0.145, Pgi=(-2.573, -3.977) 


(acceptable — needs further investigation) 


Scenario D 


85.0 


62.6 


1.9 


20.3 


0.04 


f)o=1.569, Pd=(-1.746, -1.776), 


relative magnifications are consistent 














60=0.43, 7s= (0.10, -41.5), 


with the data presented here 














6gi=0.146, Pgi=(-2.573, -3.979) 


(acceptable — needs further investigation) 


position of 


41.0 


21.3 


3.1 


16.3 


0.2 


bo=1.551, Po=(- 1.799, -1.779), 


an offset of ~ 60 mas from the optical 


Gal D free 












eo=0.43, 7s=(0.11, -41.8), 


position is significantly large 


- Scenario D 












f,gi=0.144, Pgi=(-2.581, -3.959) 


(Not acceptable) 


SIE+SlS+SIS+shear 
Scenario A 


79.5 


59.0 


1.9 


18.0 


0.7 


bo=l-551, Po=(- 1-799, -1.779), 


Adding G2 does not improve 



60 =0.42, 7s= (0.09, -40.5), the model fit to the data 

6gi=0.146, Pgi=(-2.574, -3.964) (Not acceptable) 
bG2=0.075, Pg2=(-5.749, 1.767) 



Here, Ai and Bi denote vectors in images A and B respec- 
tively, and are counterparts of each other. Thi s formulation can be 
applied to real gravitational lens systems (e.g. jGarrett et aljl994bl : 
Ijin et al.ll2003l) provided the lensed images show a rich structure of 
non-coUinear features. Note that at least two vectors are needed to 
constrain the four independent elements in the relative magnifica- 
tion matrix. 

Consider a set of three components in A (Al, A2, A3) and 
B (Bl, B2, B3). Let A12 (B12) and ^32 (B32) denote the vectors 
originating from component A2 (B2) and ending at Al (Bl) and 
A3 (B3) respectively. The positions of the components of images 
A and B are taken from Tables [5] after applying the appropriate 
shift of origin. The uncertainties are chosen as 0. 1 mas for all po- 
sitions throughout for simplicity. The determinant of the relative 
magnification matrix is found to be —1.09 ± 0.35 and is in good 
agreement with the observed flux density ratios of the images. The 
error-bars were obtained by drawing 1000 realizations of the image 
positions which follow normal distributions centred at the original 
image positions and with a scatter equal to the uncertainty in the 
positions. The negative value of the determinant again confirms the 
opposite parity of image B with respect to image A. 

We further use the results of numerical simulations to deter- 
mine the number of massive subhaloes near images A and B that 
would produce observable astrometric distortions. We use the sub- 
halo mass function from Eq. [T]to find the number of subhaloes in 
the mass range lO'^ Mq to 10^" M0. The number of subhaloes 
at the projected separation of the images, thus, expected within an 
angular region of lOOx 100 mas^ is ~ 10~^. The number of sub- 
haloes within an annulus of width 100 mas at the projected separa- 
tion of the images is ~ 1. 



7 CONCLUSIONS 

Multi-frequency high-resolution radio observations of the gravita- 
tional lens MG 2016+112 were conducted to carry out a spectral 
analysis and to find a mass model for the complex structure in the 
lensed images. Radio maps made with simultaneous MERLIN and 
global VLBI observations at both 1.7 and 5 GHz were presented. 
Subsequently, HSA observations at 8.4 GHz were undertaken to 
carry out a spectral study of the components at high resolution. In 
addition to the two previously known components in images A and 
B, three new components were detected in the observations pre- 
sented here. The observations with the HSA proved crucial in the 
confirmation of the new components. A total of five components 
are now found in images A and B. No more new components are 
detected above 33 /iJy within a region of size 0.21x0.21 arcsec^ 
centred at images A and B from the HSA imaging at 8.4 GHz. A 
5a upper limit was placed on the peak surface brightness of an odd 
image in the vicinity of the lens D, or radio emission from D, of 
0.18 mjy beam~^ at 8.4 GHz. 

The overall radio spectra and the flux densities of the compo- 
nents in A and B were found to be similar. The flux density ratio 
of images A and B were also found to be consistent with the de- 
terminant of the relative magnification matrix. Therefore, there is 
no significant substructure or any other effects that might affect the 
flux densities of the images. In region C, the morphology and spec- 
tra of CI I-C21 and C12-C22 were found to be similar, as expected 
for lensed images. Furthermore, the measured flux-density of the 
C2 pair is ~1.2 times the flux-density of the CI pair at all frequen- 
cies, which could be due to the proxi mity of the sa tellite galaxy Gl 
to the C2 pair with a positive parity ( lKeetonll2003b . The identifica- 
tion of components in region C with those in image A (or B), on the 
basis of their spectra, cannot be done because the highly magnified 
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components of region C correspond to extremely small regions in 
either the detected (4, 1, 3, 2 or 5) or undetected components of 
images A and B. 

Several mass models with more than one mass component in a 
single lens plane were investigated for four scenarios. In these sce- 
narios, the components of region C were constrained as the lensed 
counterparts of different parts of the components of images A and 
B, and the consequences of doing so were assessed. The mass mod- 
els tested here with scenarios A and B predicted relative magnifi- 
cations of the images that were inconsistent with the observations. 
The predictions of scenarios C and D were consistent with the ob- 
servations presented he re. The predictions of the mass models of 
iKoopmans et al.l j2002l) are not consistent with the new observa- 
tions because of the detection of component 5. A SIE+SIS+shear 
model with the satellite galaxy Gl (SIS) found at the same redshift 
as the lensing galaxy D (SIE) improved the fit (reduced = 3.5) 
to the astrometric constraints significantly, a s stated previously 
bv lKochanek. Schneider & WambsganssI j2004) . However, a model 
with more complexity, for example, ellipticity in the satellite galaxy 
which is currently observationally ill-constrained but can cause lo- 
cal azimuthal distortions near region C, is perhaps required. 

We have compared the subhalo abundance and mass fraction 
in MG 2016+112 with lensing and simulations. About 1 to 2 sub- 
haloes are expected in an annulus around the images from the 
CDM scenario and we detect one (luminous) subhalo centred at 
Gl. Prima facie, this appears like an agreement which can be ei- 
ther confirmed or refuted with further deeper observations of the 
environment of this lens system. The subhalo mass fraction, as- 
suming an isothermal profile, indicates a much higher mass in the 
su bhalo centred at G l than that expected from the simulated halo 
of iMao et af] j2004l) . Considering the assumptions and caveats in- 
volved here, this comparison should only be taken as suggestive 
rather than conclusive. 
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APPENDIX A: DEFINITIONS AND EQUATIONS 

The critical radius of an isothermal ellipsoidal lens is related to the 
velocity dispersion {a) through the following equation 



■ fcsis = \ , , , 47r — — - (Al) 



which takes into account the ellipticity through the axis ratio q. The 
mass of the lens within the Einstein radius can then be calculated 
from 

M«.s:s)^%f^. (A2) 

Here, Dds, and Ds refer to the angular diameter distances be- 
tween the lens and the source, the observer and the lens, and the 
observer and the source, respectively. 

The definition of virial radius ryjr and the mass Mvir of the halo 
within the virial radius jMo, Mao & WhitJ 19981) is given by 

Vc 

(A3) 



10 H[z) 

= 10^ • 

where Vc = \/2a is the circular velocity, G is the gravitational 
constant and H{z) is the Hubble parameter. 

The mass of satellite Gl within 0.025 rvir, the mass of lens galaxy 
D within an annulus and the satellite mass fraction are calculated 
using the following equations, 

«»■ - L '"'^ 

_ Mgi _ (^Gi 0.025 ryj, 
^^"^ " Md ' R1-R2 ^ ' 

Here, Rl and R2 correspond to the projected radial distances 
from the center of the potential which bound the annulus. 



